Quantitative test of thermal field theory for Bose-Einstein condensates II 
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We have recently derived a gapless theory of the linear response of a Bose-condensed gas to 
external perturbations at finite temperature and used it to explain quantitatively the measurements 
of condensate excitations and decay rates made at JILA [D. S. Jin et ai, Phys. Rev. Lett. 78, 764 
(1997)]. The theory describes the dynamic coupling between the condensate and non-condensate via 
a full quasiparticle description of the time-dependent normal and anomalous averages and includes all 
Beliaev and Landau processes. In this paper we provide a full discussion of the numerical calculations 
and a detailed analysis of the theoretical results in the context of the JILA experiment. We provide 
unambiguous proof that the dipole modes are obtained accurately within our calculations and present 
quantitative results for the relative phase of the oscillations of the condensed and uncondensed atom 
clouds. One of the main difficulties in the implementation of the theory is obtaining results which are 
not sensitive to basis cutoff effects and we have therefore developed a novel asymmetric summation 
method which solves this problem and dramatically improves the numerical convergence. This new 
technique should make the implementation of the theory and its possible future extensions feasible 
for a wide range of condensate populations and trap geometries. 

PACS numbers: 03.75.Kk, 67.40.Db, 05.30.Jp 
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I. INTRODUCTION 



Measurements of the excitation spectrum of Bose- 
Einstein condensates (BECs) provide a unique opportu- 
nity to compare quantitative predictions of finite tem- 
perature quantum field theories (QFTs) with experimen- 
tal data in a regime where interactions, finite temper- 
ature and particle statistics are all important. At low- 
temperature the experimental results can be understood 
using only the Gross-Pitaevskii equation (GPE) for 
the condensate, but the situation is less satisfactory at 
finite temperature where a description of the coupled dy- 
namics of condensed and uncondensed atoms is required. 
There are also more fundamental difficulties with de- 
veloping a consistent finite temperature theory, such as 
dealing correctly with the ultra-violet and infrared diver- 
gences which can appear and obtaining a gapless spec- 
trum, as required by the Hugenholtz-Pines theorem 0. 
We have recently developed a gapless theory which ad- 
dresses all these issues and provides a general method for 
calculating the response of a BEC to external perturba- 
tions at finite temperature [3|- We have demonstrated 
the validity of this approach by applying it successfully 
to the pioneering measurements of condensate excitations 
made at the Joint Institute for Laboratory Astrophysics 
(JILA) in 1997 0, 0. In particular, we were able to 
explain the sudden upwards shift in the resonance fre- 
quency of the low-lying m — mode as the tempera- 
ture increased towards the critical temperature for BEC 
formation. In this paper we describe our numerical cal- 
culation and provide a detailed analysis of the results 
and their implications in the context of this experiment. 
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The theoretical approach and numerical techniques are 
general, however, and can be applied directly to a wide 
variety of recent experiments on BECs [E 0> HI @ • 

The early experiments at JILA in 1996 and 1997 |ilicj 
measured the energies and decay rates of low-lying con- 
densate excitations with axial angular momentum quan- 
tum numbers m = 2 and m — as a function of 
the condensate population and temperature. The low- 
temperature measurements are in good agre ement with 
the usual Bogoliubov quasiparticle theory |loL fTTL and 
with calculations based on the Hartree-Fock-Bogoliubov 
(HFB) formalism but the finite-temperature results 
have proved much harder to explain. In the experiment 
of the m = 2 mode was observed to shift strongly 
downwards with temperature, while the m — mode 
underwent a sharp increase in energy towards the non- 
interacting limit. The results for the m = 2 mode could 
be explained using a gapless extension of the HFB ap- 
proach (GHFB) which includes the so-called anomalous 
(pair) average of two condensate atoms 0, flfij . and 
also by the dielectric formalism [lfi| . However, both ap- 
proaches were unable to explain the upward shift of the 
772 = mode, and an analytical calculation also predicted 
downward shifts for both modes ^3 • 

A possible explanation for the behaviour of the m = 
mode was given by Bijlsma, Al Khawaja and Stoof 
[Tsl llflj in terms of a crossover from out-of-phase to 
in-phase oscillations of the condensed and un-condensed 
atoms at high temperature. Subsequently, Jackson and 
Zaremba [2(J] obtained good agreement with the JILA 
results for both modes using the Zarcmba-Nikuni-Griffin 
formalism 0, 122I l2^ | which involves coupling a general- 
ized GPE for the condensate to a a semi-classical Boltz- 
mann equation for the non-condensate. Despite its suc- 
cesses, however, this approach neglects the anomalous 
average and all Beliaev processes. Although Beliaev pro- 
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cesses are expected to be swamped by Landau processes 
at high temperature, they have nonetheless been directly 
observed in a number of recent experiments EH^EHEa], 
while the good agreement with the JILA results for the 
m = 2 mode within the GHFB theory 0] suggests that 
the anomalous average can also be important. 

In two recent papers, we have developed a systematic 
perturbative extension of the Bogoliubov theory which 
includes these effects and explains the JILA results for 
both the m = 2 and m = modes |3|, la] • The formal- 
ism adapts the linear response treatment of Giorgini 0] 
and provides a time-dependent extension of an earlier 
second-order perturbative calculation j^. The theory 
is gapless and includes the dynamic coupling between 
the condensate and non-condensate, all relevant Beliaev 
and Landau processes and the anomalous average. It is 
also consistent with the generalized Kohn theorem for the 
dipole modes 0]. The theory is valid in the collision- 
less limit of well-defined quasiparticles, which requires 
(fc B T/n C/o)(noa3) 1 /2 < i whe re n Q is the condensate 
density, a s is the s-wave scattering length, fee is Boltz- 
mann's constant and Uq = 4irh 2 a s /m where m is the 
atomic mass ^j, E2 ■ For the JILA experiment £| this 
parameter does not exceed 0.03 at the trap centre for the 
highest temperature we consider. 

An important feature of our theoretical approach is 
that it explicitly includes the effect of the external per- 
turbation on the non-condensate dynamics. At high tem- 
peratures and for perturbations which are peaked on the 
edge of the atomic cloud, the non-condensate has a large 
response arising from single-particle resonances at inte- 
ger values of the trap frequencies. If the quasiparticle 
mode is located close to such a resonance, it is possible 
to excite the condensate via the thermal cloud as inter- 
mediary and the strength of this excitation can exceed 
the direct effect of the external perturbation. This turns 
out to be the explanation for the anomalous behaviour of 
the m = mode in the JILA experiment |a] and it is also 
necessary to include this process for a correct description 
of the dipole oscillations of the system, as we show later 
in this paper. 



A. Outline of the paper 

This paper is organised as follows. In Sec. [D] we sum- 
marize the relevant results of the theory developed in 
Ref. |3 and in Sec. II I II and its subsections we discuss 
our numerical methods. We focus particularly on the 
important issue of convergence and in Sec. IIII El we de- 
scribe a new asymmetric summation technique we have 
recently developed which greatly improves the conver- 
gence properties of the calculation. In Sec. lIVI we present 
our numerical results for the parameters of the JILA ex- 
periment. We focus first on the results for the m = 2 
and m = modes, describing the effect of some of the 
Landau and Beliaev processes which can take place and 
showing the importance of including the effect of the per- 



turbation on the thermal cloud. In Sec. IIV CI we anal- 
yse the experimental results for the condensate popu- 
lation and temperature to determine the correct input 
parameters for the theory and in Sec. IIV Dl we provide 
results for the dipole modes and show that these are ob- 
tained to high accuracy in our calculation. Finally, in 
Sec. IIV El we present results for the relative phase of the 
condensate-non-condensate oscillations and demonstrate 
that the cross-over from out-of-phase to in-phase oscilla- 
tions predicted for the m = mode at hig h temperature 
by Al Khawaja, Bijlsma and Stoof 0,0] is reproduced 
in our theory. Two appendices contain details of the cal- 
culation and the definitions of some of the quantities we 
require. 



II. SUMMARY OF SECOND ORDER 
BELIAEV-POPOV THEORY 

In this section we briefly summarize the theory that 
we use in this paper, which we will refer to as the second 
order Beliaev-Popov (SOBP) theory. The full derivation 
of the expressions we quote along with a detailed discus- 
sion of their physical significance is given in Ref. 0, to 
which we frequently refer. 



A. Equations of motion and Bogoliubov 
quasiparticles 

The SOBP theory describes the evolution of a conden- 
sate in the presence of non-condensed atoms using a gen- 
eralized GPE containing various non-condensate mean- 
fields. These are calculated by modelling the thermal 
cloud as a non-interacting gas of quasiparticles evolving 
in the time-dependent mean-field of the condensate ac- 
cording to an appropriate set of Bogoliubov-de Genncs 
(BdG) equations. Together the GPE and BdG equations 
form a consistent description of the coupled dynamics of 
the system which contains all the leading order correc- 
tions to the usual Bogoliubov theory in the relevant small 
parameter, defined above. 

The generalized GPE for the condensate wavefunction 
<l>(r,i) (normalized to one) has the form 



i»|*(r,t) 



H sp (r) + P(v,t) - X(t)\ $(r,t) 
[N (t)UoMr,t)\ 2 + 2U n(r,t)] $(r,i) 
lWHr,i)$*(r,i)-/(r,t), (1) 



where H sp {r) = — ft 2 V 2 /2m + I4ra P (r) is the single- 
particle Hamiltonian containing the static trap potential 
(if any), P(r,t) is a time-dependent external perturba- 
tion, A(i) is a scalar which controls the global phase of 
$(r, t), and No(t) is the condensate population 29]. Par- 
ticle interactions are described using a contact potential 
with interaction strength Uq- The use of this contact po- 
tential necessitates an ultra-violet (UV) renormalization 
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of the theory, which is achieved by a suitable modification 
of the pair average m R (r, t) (see below). 

The coupling of condensed and non-condensed atoms is 
described by the terms involving the non-condensate den- 
sity h(r,t), the renormalized anomalous (pair) average 
m R (r, t), and the function f(r,i) which arises from the 
fact that the wavefunctions of the non-condensed atoms 
are orthogonal to the condensate. These quantities are 
constructed from a complete set of time-dependent quasi- 
particle wavefunctions ttj(r, t) and Wj(r, t) according to 



n(r, t) 



]T Mr, i) | 2 TV, + \v t (r,t)\ 2 (N % + 1), (2) 



m R (r, t) = £ Ui (r, t)v*(r, t)(2N t + 1) + ^^$ 2 (r, t), 



?7n 



(3) 



*) = W J2 c * N Mr, t) + Ci{Ni + l)v*(r, i), (4) 



dit) = N U dr\<S>\ 2 [<f>*u l (r,t) + <5>v i (r,t)}. (5) 



The final term in the expression for ?ri R (r, t) is the UV- 
renormalization and the quantity AUq is the second-order 
approximation to the interaction strength Uq as calcu- 
lated from the Lippmann-Schwinger equation 



AU = U Q 2 



d 3 k 



1 



(2tt) 3 2(H 2 k 2 /2m) 



(6) 



To obtain a closed set of equations we also require the 
equation of motion for the quasiparticle wavefunctions. 
These evolve according to the BdG equations 
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L(r, t) = H sp + P(r, t) + N U [|$| 2 + Q|$| 2 Q 
M(r,t) = N U Q<Z> 2 Q*, 



(7) 

(8) 
(9) 



where the projector Q = 1 — |$) ($| ensures orthogonality 
of $(r, t) and {ui(r 7 t),v* (r,t)}. The quasiparticle pop- 
ulations {iVj} are independent of time and given by the 
Bose-Einstein distribution Ni = \/[e l3t ' €i ~ 5 ^ — 1), where 
d is the Bogoliubov energy (see below), j3 = 1//cbT and 
Sfi is the (small) difference between the condensate en- 
ergy and the chemical potential. The condensate pop- 
ulation is defined implicitly using the constraint on the 
total particle number N by No(t) — N — Jdr n(r, t). Most 
quantities in the theory depend on temperature via their 
dependence on the quasiparticle populations. 

The above equations are obtained using the number- 
conserving formalism of Gardiner and Castin and Dum, 
modified for finite temperature calculations 0, l30l l3ll 
l32l |. The terms /(r, t) and Q which arise from orthogo- 
nality of the condensate and non-condensate are a feature 
of this approach and do not appear in symmetry-breaking 
theories. We find numerically that they can give a sig- 
nificant contribution to the energy shifts (see Sec. If 11 F|l . 



B. Equilibrium solutions 

In equilibrium, P(r, t) = and Eq. (JTJ has a time- 
independent solution $(r,t) = <i>(r) which satisfies 

H sp (r) — A + AWo|*(r)| 2 + 2C/ n(r)] $(r) (10) 
+C/ m R (r)$*(r)-/(r) = 0, 

where A is the condensate eigenvalue and n(r), m R (r) and 
/(r) are equilibrium non-condensate mean-fields calcu- 
lated from Eqs. using the static quasiparticle basis 
defined below. Setting these quantities to zero gives the 
usual time-independent GPE with wavefunction < E'o( r ) 
and energy Aq 



ff sp (r)-A + AWo|$o(r)| 2 $o(r) = 0. (II) 



We solve Eq. (|1(J|> by linearizing the change in energy and 
shape of the condensate relative to this solution. 

We obtain time-independent BdG equations by writ- 
ing $(r,i) -> $ (r) and P(r,t) = in Eqs. 0-@ and 
looking for solutions of the form 



Ui(r, t) 
«i(r, t) 



Ui{r) 
Vi(r) 



(12) 



The equilibrium Bogoliubov quasiparticle wavefunctions 
and energies are therefore found by solving 



La M 
-M * -L* 



(13) 



where Lq and Mq are defined by making the appropriate 
substitutions in Eqs. l(5 )l -l(5 |l . An important consequence 
of the projector Q in Lq and Mo is the existence of two 
zero-energy solutions proportional to $o ( r ) which means 
that these quasiparticle wavefunctions form a complete 
basis set. In particular, we obtain a set of solutions 



(14) 



with energies ej,— £i,0, and norms +1,-1, +1,-1 re- 
spectively, where the norms are defined by Jdr\ui\ 2 — 
\vi\ 2 . We use the notation i = 0+ to denote the positive- 
norm, zero-energy mode, which we will also refer to as 
'the zero mode'. The ease with which the two zero energy 
solutions are obtained and used in the theory is an im- 
portant advantage of the number-conserving approach. 



C. Linear response theory 

We can use the above equations to find the linear 
response of a Bose-condensed gas to an external per- 
turbation. We consider the situation where a conden- 
sate has been formed at low temperature and has set- 
tled into the ground state of Eq. ljTU)) 33|. When the 
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external perturbation P(r, t) is applied, the system re- 
sponds with a time-dependent oscillation of all mean 
fields around their static values, $(r, t) — $(r) + <5$(r, t), 
n(r, t) = n(r) + Sh(r, t) etc. Substituting these expres- 
sions into Eq. Q and linearizing leads to the equation 
of motion for the condensate fluctuation <5$(r,<). This 
equation can be solved by combining it with its com- 
plex conjugate, Fourier transforming and expanding in 
the equilibrium quasiparticle basis 



(5$*(r, —a/ 



Ui(r) 
Vi(r) 



(15) 



The expansion coefficients bi{uj) (response amplitudes) 
are directly related to the condensate density fluctuations 
Sn c = S(N \^\ 2 ), which are the experimentally relevant 
quantities 



6n c {r,Lj)=6N (w)\$ \ 



(16) 

where SNq(u>) = — Jdr Sn(r, oj) describes any fluctuations 
in the condensate population. This quantity is only non- 
zero if the perturbation has a contribution with the same 
symmetry as the condensate (as for the m = mode in 
the JILA experiment for example) and we have found 
numerically that it generally has rather little effect on 
the observed density fluctuations. 

In general, the expansion coefficients bi{u) can be 
found from the solution of a linear matrix problem. The 
result is particularly simple, however, if a single positive- 
norm mode dominates the expansion, which is usually 
the situation for experiments designed to study excita- 
tions. In the case that mode 'p' is excited, the solution 
has the form 



(17) 



where 7 measures the experimental resolution, the exci- 
tation matrix element P p o(lu) is defined by 



JdrP(r,co) [«;* + «;*o] 



(18) 



and the response function (resolvent) T p {uj) takes differ- 
ent forms depending on the level of approximation. In the 
simplest case where the coupling to the non-condensate 
is completely neglected, F p {uj) is a lorentzian centred on 
the frequency of the corresponding Bogoliubov mode 



T p {u) 



two — 



(19) 



This function diverges at the resonance frequency uj p — 
e p /H because there is no intrinsic damping in the the- 
ory at this level of approximation. The inclusion of the 
resolution parameter 7 in Eq. i|17[) via the standard sub- 
stitution u) — > lu + 27 ensures we obtain finite quantities 
(essential for numerical work) and can be justified from 
the finite experimental observation time [j| . Typically, 7 



is of order a few hundredths of an appropriate trapping 
frequency and our results do not depend strongly on its 
precise value within the experimentally relevant range. 

If we now include the dynamic coupling between the 
condensate and non-condensate, the response function 
becomes 



T p (uj) = U p (lu), 



AP( S) M + AP^>) 



pO 



(D) 



pO 



-Ppo(w) 



G P {u) 



G P {u) ■ 
E p (l>) 



huj — E p uj) ' 



T, p (oj). 



(20) 
(21) 

(22) 

(23) 
(24) 



Here E p (w) is a frequency-dependent self-energy, while 

the quantities AP^\uj) and AP^\oj) depend on the 
functional form of the external perturbation and describe 
changes in the excitation amplitude P p $(uj) resulting from 
the coupling of the condensate to the thermal cloud. 

The various response functions introduced above de- 
scribe different dynamical processes occurring in the gas. 
In general both the condensed and uncondensed atoms 
can be excited via two distinct mechanisms; either di- 
rectly by the external perturbation or indirectly by fluc- 
tuations in the other mean-fields. For the case of the 
condensate, these two mechanisms are described by the 
response functions G p (lo) and G p (uj) respectively, while 
the total response including both processes is described 
by fcp{u>). 

The separate response functions Gp(u) and G P {u) are 
introduced partly to facilitate interpretation of the the- 
ory and partly because they may dominate in certain sit- 
uations. In general, the relatively large population and 
density of the condensate means that the dominant pro- 
cess at low temperature is for the perturbation to ex- 
cite the condensate which subsequently drives the non- 
condensed atoms via their mutual coupling. This case 
is described by Gp(w) alone. The alternative mechanism 
where the perturbation excites the non-condensate first 
and this acts on the condensate in a second step is de- 
scribed by Gp{w)- Under certain circumstances, (and es- 
pecially at finite temperatures) this second process can 
become dominant and its inclusion is crucial to explain 
the JILA experiment and for an accurate description of 
the dipole modes, as we show in Sec. HVDl 0, 0. It 
is also possible to enhance the effect of one or other of 
these mechanisms by a suitable choice of the perturbing 
potential P(r,t). If this is localized on the condensate 
then the response will be dominated by Q p (ui), while if 
it mainly acts in the wings of the non-condensate then 
G P (u) is more appropriate. This later case has been used 
at MIT to investigate second sound oscillations 0. We 
stress, however, that the full response is given by 1Z p (uj) 
and includes both mechanisms. 

The self-energy in Eq. Q24[l contains two distinct types 
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of contributions, static (S) and dynamic (-D), correspond- 
ing to the different roles of the thermal cloud 

X p {w) = AEjj?+AEg\w). (25) 

The frequency-independent static term AE PP comes 
from interactions between a condensate fluctuation and 
the equilibrium non-condensate mean-fields, while the 
dynamic term AE P p^ (u>) describes the driving of the non- 
condensed atoms by the condensate and their subsequent 
back action. 

The detailed definition of these quantities is given in 
3] and in Appendix 1X1 Briefly, the static term can be 
written as a sum of contributions of the form 

AE^ = AE 4 ( P ) + AE X ( P ) + AE shapc (p) + AE { /\p). 

(26) 

Here AE^ip) and AE^ (p) arise from the explicit inter- 
actions between a condensate fluctuation and the equi- 
librium mean-fields h(r), m R (r) and /(r), while AE\(p) 
and A^shapcCp) come from the effect these mean- fields 
have on the energy and shape of the equilibrium conden- 
sate [i.e. the difference between the solutions to Eqs. <|f Of) 
and All static terms can be calculated straight- 

forwardly from integrals involving the condensate wave- 
function and the equilibrium non-condensate mean-fields. 
The dynamic terms on the other hand involve a double 
summation over the quasiparticle basis states of the form 

A4?>,Aiff>~£ f{Ni)M ™ . (27) 

ij J 

Here f(Ni) is a simple function of the quasiparticle pop- 
ulations Ni, Mpij is a suitable matrix element and w%j 
is a resonance of the non-condensate corresponding to 
a Beliaev or Landau process (wy = ±(e, + ej) and 
u>ij = ±(ej — Ej) respectively). The matrix elements Mpij 
are the product of two factors each of which is defined 
in terms of integrals of the equilibrium condensate wave- 
function < &o( r ) and the quasiparticle wavefunctions for 
modes p, i, and j. 



III. NUMERICAL METHOD 

We have used the formalism outlined above to calcu- 
late the density response of a trapped condensate as a 
function of temperature. In this section we describe the 
numerical techniques required and present some illustra- 
tive results. The calculation is difficult because some of 
the intermediate terms we require are much larger than 
the final self-energy so there is substantial cancellation 
between them. This puts significant demands on numer- 
ical accuracy and convergence which we have solved using 
a variety of techniques, described briefly below. Further 
details of the methods used can be found in Ref. |34j . 



We consider the case that the trapping potential is an 
axisymmetric, anisotropic harmonic potential of the form 

H rap (r) = im(^r 2 +a;^ 2 ), (28) 

where r 2 — x 2 + y 2 is the square of the radial coordi- 
nate. In the present paper we use simulation parame- 
ters appropriate to the 1997 JILA experiment [4]. The 
trap frequencies are therefore v r — uj r /(2ir) = I29Hz and 
v z = uo z / (2tt) = 365Hz, the scattering length for the Rb 87 
atoms is a s = II0a o where clq is the Bohr radius and the 
resolution parameter 7 = 0.36w r These parame- 

ters are fixed for all the results presented in this paper. 
In addition we take the condensate population No to be 
6000 unless specifically stated otherwise. 

For a given atomic species and trap geometry, the vari- 
able input parameters of the theory are the condensate 
population Nq and the absolute temperature T. For a 
given value of Nq, the numerical calculation of the con- 
densate density response can be broken down into the 
following steps: 

1 . Solve the time- independent GPE of Eq. ifTTl) to ob- 
tain the static condensate wavefunction <&o( r ) and 
eigenvalue Ao- 

2. Solve the equilibrium BdG equations of Eq. I|13l) to 
obtain the static quasiparticle wavefunctions Uj(r) 
and Vi(r) and energies e, for all states up to a nu- 
merical cutoff energy E cut . These solutions are 
saved to disk along with the condensate wavefunc- 
tion and eigenvalue. 

3. For each temperature T, construct the equilib- 
rium non-condensate density n(r) and renormalized 
anomalous average m R (r) on a spatial grid by a 
direct summation over the quasiparticle states us- 
ing the time-independent form of Eqs. J3J and J3J), 
supplemented by a semi-classical approximation for 
states above the numerical cutoff (see Appendix^ . 
From these, the equilibrium solutions to the gener- 
alized GPE can be found straightforwardly, as can 
all the static shifts defined in Eq. (|26[) . These calcu- 
lations simply involve integrations over the spatial 
grids defining the equilibrium mean-fields. 

4. For each mode 'p', temperature T and resolution 7 
the dynamic terms AE p D \uj) and AP^\lo) (and 
8Nq{uj) if required) are calculated on a frequency 
grid centred on the Bogoliubov energy e p . This 
requires looping over the double quasiparticle sum- 
mation, calculating the non-zero matrix elements 
for the relevant Landau and Beliaev processes and 
combining these with the appropriate energy de- 
nominator. The magnitude of the task is greatly 
reduced by selection rules which mean that only a 
relatively small fraction of the terms in the sum are 
non-zero. For AP„q (w), this step requires a choice 
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for the functional form of the external perturbation 
P(r, t), and this is taken to mimic the experiment 
as closely as possible [see Eq. (j2HJ)]. In some cases a 
semi-classical approximation is also used to include 
the effect of high energy states above the numerical 
cutoff. 

5. The various contributions are combined to give the 
self-energy Ti p (u> + 17) and the response functions 
Q p (uj + ry), G p (u + iy) and 1Z p (uj + ij). Energies 
and decay rates are extracted from these quantities 
either by reading off the value of the self-energy at 
the poles or by fitting appropriate functions to the 
resolvents in either the time or frequency domains. 

6. To compare with the experimental data, our results 
must be plotted against the reduced temperature 
t = T/T® rather than the absolute temperature T. 
Here T® is the critical temperature of an ideal gas 
of N atoms in a harmonic trap, given by 



N 

W) 



1/3 



(29) 



where u> = (o^u^) 1 / 3 is the geometric mean trap 
frequency, C( x ) is the Riemann zeta function, and 
C(3) « 1.202. The conversion requires calculat- 
ing the total atom number N for each temper- 
ature, which is achieved using the relation N — 
No + Jdr n(r). This gives N, Tj? and t as a function 
of temperature in nanokelvin for a given fixed value 
of the condensate population Nq. 

Typically, we solve the BdG equations with a very 
large energy cutoff of E cut ~ 14Qhw r which corresponds 
to more than 100, 000 quasiparticle modes. This takes 
about 5 hours of computation on a desktop PC with a 
2.4GHz Pentium 4 CPU and 1Gb RAM and consumes 
just over 2Gb of storage on disk. The subsequent cal- 
culation of all the relevant response functions for 20 dif- 
ferent temperatures and a frequency grid of 512 points 
takes roughly an additional 5 hours of computation for 
each excitation under study. In the present work these 
are restricted to the lowest energy m — 2 and m = 
modes and the dipolc oscillations. 

In the following sections, we provide further details on 
some of the above steps along with illustrative numerical 
results. 



A. GPE solution 

The GPE of Eq. (|llfl is solved by expanding the wave 
function < I > o( r ) hi an appropriate basis and solving the 
resulting set of coupled nonlinear equations for the ba- 
sis coefficients. The basis states are chosen to be prod- 
ucts of the harmonic oscillator (HO) states for the radial 
and axial directions, while the solution of the nonlinear 
equations is achieved using the MATLAB optimisation 



routine FS0LVE, with an initial guess provided by the 
Thomas-Fermi solution. This is an efficient method of 
solving the GPE as it is relatively stable initially while 
also giving rapid (quadratic) convergence near the end of 
the solution cycle. Various integrals of products of four 
functions (such as the projection of l^oP^o onto a basis 
state for example) are required to set up the nonlinear 
problem. These are calculated on a Gaussian quadrature 
grid chosen to allow exact integration of such products. 



B. BdG solution 

The equilibrium BdG equations given in Eq. (1131) are 
solved using the procedure described in Ref. 35] . The 
symmetry of the trap potential means that the equations 
decouple into subspaces with definite values for the z- 
component of angular momentum m and z-parity p = 
(even) or 1 (odd). Within each subspace, the solutions 
are assigned another quantum number n which orders 
the energy within that subspace. The subscript 'i' on 
the quasiparticle wavefunctions Ui (r) and Vi (r) therefore 
stands for the triplet of quantum numbers (mi,Pi,rii). 
The quasiparticle energies and the radial and axial parts 
of the wavefunctions only depend on the modulus of m; 
so we only need to solve the BdG equations in subspaces 
with mi > 0. The dependence on the angular coordi- 
nate 4> has the usual complex exponential form and the 
solutions can therefore be written as 



Ui(r) = W|roil>P«>»»( r > z )' 



irrii4> 



2tt 
/2¥ ' 



(30) 
(31) 



Within a given subspace of m and z-parity, the BdG 
equations are solved in two stages. In the first stage, we 
obtain a single-particle basis orthogonal to the conden- 
sate and with the appropriate symmetries by solving the 
equation 



# sp (r) + AW |$oM)| 2 CiW 



'Ci(r). (32) 



We will refer to the solutions of this equation with 
e' 6 ^' ^ A as the GPE basis. Equation is solved 
by expanding the wave functions {Q} using the under- 
lying HO basis set. The BdG equations are solved by 
rewriting them as decoupled equations for Ui ± Vi and 
expanding in the GPE basis which reduces them to a 
standard matrix eigenvalue problem. Finally, the solu- 
tions are converted from the GPE basis to the HO basis 
and stored. The GPE basis wave functions are intro- 
duced for two reasons; firstly, they are orthogonal to the 
condensate and therefore the treatment of the orthogo- 
nal projector is trivial and secondly the BdG equations 
take a particularly simple form when written in the GPE 
basis 351. 
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We solve for all quasiparticle modes with energies 
less than a numerical cutoff E cut determined by set- 
ting a maximum value for the angular momentum quan- 
tum number m, which we denote by m max . E cut is 
then taken as the lowest energy state in the subspace 
with m = m max and even z-parity. In our largest cal- 
culations we take m max = 150 which corresponds to 
E cut « 140fiti> r . Such a large basis set is necessary to 
ensure numerical convergence of our final results (see 
Sec. 1111 E[) and requires that care be taken to maintain 
numerical accuracy at all stages of the calculation. 

The solution of the BdG equations requires the calcu- 
lation of numerous integrals involving products of either 
two or four functions. Integrals of two functions can be 
done exactly using the scalar product of the HO expan- 
sion coefficients. Integrals of four functions are calculated 
on a Gaussian quadrature grid, as in the solution of the 
GPE. The grid points and weights are chosen to allow 
exact integration of any product of four functions con- 
structed from the underlying HO basis. However, since 
all the integrals required contain the condensate density, 
we discard all points which lie outside the region where 
this has fallen to 10~ 10 of its peak value. This has es- 
sentially no effect on the accuracy of the integration but 
keeps the numerical grid to a reasonable size for calcula- 
tions with a large basis. 

C. Matrix elements for dynamic terms 

The dynamic terms require the calculation of numer- 
ous matrix elements and these are produced 'on-the-fly' 
by loading the stored quasiparticle wavefunctions in the 
HO representation, constructing them on a suitable spa- 
tial grid and summing over the points. The number of 
integrals which must be evaluated is kept manageable us- 
ing angular momentum and z-parity selection rules which 
mean that only a small subset of the matrix elements 
are non-zero. For the dynamic self-energy AE p D ^ (u>), 
the integrals involve products of four wavefunctions with 
two of low energy (the condensate wave function and the 
quasiparticle mode 'p' under study) and two with ener- 
gies up to the numerical cutoff E cu t . These integrals can 
therefore be done essentially exactly using the truncated 
Gaussian quadrature grid used in the solution of the BdG 
equations. 

Unfortunately, this grid can not be used to calcu- 
late the matrix elements required for AP^'(w) which 
involve the product of the perturbation P(r, u>) in the 
space-frequency domain and two quasiparticle wavefunc- 
tions. The reason is that P(r, oj) is generally not local- 
ized in the centre of the condensate but extends out to 
the wings of the trapped cloud. However, the perturba- 
tion typically involves only very low order polynomials 
so these integrals can instead be evaluated by construct- 
ing a new Gaussian quadrature grid which can integrate 
bilinear products of quasiparticle functions exactly and 
then adding the few extra points required to deal with 



the perturbation. In this way all integrals required in the 
theory are calculated exactly and the accuracy of the nu- 
merics is limited purely by the size of the basis employed 
(i.e. by the value of E cut and, for a given E cut , by the 
size of the HO basis used to construct the solutions). 

The functional form of the perturbation is chosen to 
mimic the experiment as closely as possible. In general, 
we take 

P(r, t) = -P p (r, z) cos{uj d t - m p (/))Q{t)&(t d - t), (33) 

where ojd is the central drive frequency, O(i) is the unit 
step function, and the perturbation is applied for < t < 
td {td ~ 14ms in the JILA experiment Q). Following the 
form of the perturbation used experimentally, we take 
P p (r, z) oc r 2 for the m = and m = 2 modes, while 
P p (r, z) oc z for the dipole mode along the z-axis and 
P p (r, z) oc r for the dipole modes in the x-y plane. 

D. Dynamic Summations 

The dynamic terms are calculated by performing the 
double sum over intermediate quasiparticle states i and j 
in Eq. (|27() . The calculation is only feasible because the 
selection rules on the matrix elements mean that only a 
small fraction of the total number of intermediate pairs 
contribute in the sum. In fact, for a given quasiparticle 
mode 'p', each value of the angular momentum rrii and 
z-parity pi for mode i corresponds to unique values for 
nij and pj, so the calculation can be broken into sub- 
space blocks. For each block, the quasiparticle energies 
and wavefunctions are loaded and the relevant matrix 
elements are calculated as described above. The sum- 
mations over the remaining quantum numbers rii and rij 
(describing the energies within each block) are then per- 
formed for each temperature T and resolution parameter 
7 and for each value of uj on a frequency grid centred 
on the energy of the mode under study. This process is 
repeated for all contributing subspace blocks and for all 
modes 'p' of interest. 

The efficiency of this calculation can be greatly in- 
creased by noting that the frequency dependence of the 
dynamic terms generally consists of a few strong features 
superposed on a smooth background. A pair of interme- 
diate states i and j is associated with an energy resonance 
for a Beliaev and Landau process when hm — ±€j ± tj [cf. 
Eq. (|27|l ] . If this resonance occurs within the frequency 
range of experimental interest, then it potentially cor- 
responds to a sharp feature which in principle we need 
to resolve. This situation therefore requires a fine grid 
spacing, although generally the sum of many such terms 
is much smoother than each individually. The fine grid is 
needed, however, if a few modes dominate the response, 
a situation which does occur in finite systems (as our re- 
sults demonstrate) and which has been termed 'temper- 
ature induced resonances' by Guilleumas and Pitaevskii 
,36]. However, the vast majority of resonances fall out- 
side the frequency range of direct interest (for example 
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there is usually at most one Beliaev resonance within this 
range for typical values). The contribution from these 
processes is smooth in frequency and can be adequately 
represented on a rather coarse grid. We therefore divide 
the intermediate states into two groups; those with res- 
onances in the range of interest are calculated on a fine 
grid while those with resonances outside this range are 
calculated on a much coarser grid and interpolated onto 
the fine grid at the end of the calculation. This scheme 
dramatically reduces the computational time at a negli- 
gible cost in accuracy. In the calculations presented here 
we used a fine grid of 512 points over a frequency range 
of 2uj r with a coarse grid of only 16 points over the same 
range. 



E. Convergence 

An important requirement of any numerical calcula- 
tion is that the results have converged sufficiently that 
meaningful conclusions can be drawn. This is particu- 
larly difficult in the present case because of residual ef- 
fects of the infrared divergence problem which plagues 
theories of the Bose gas beyond the Bogoliubov approxi- 
mation. Although the full theory is infrared finite, indi- 
vidual terms in the self-energies are not; for example, in 
the large volume homogeneous limit both the static and 
dynamic contributions diverge as 1/k as k — > and only 
their sum is finite and proportional to the small param- 
eter of the theory [H l27l |37| . 

In a finite system, this divergence is suppressed but the 
static and dynamic terms may still become large com- 
pared to other energy scales in the problem. This is 
demonstrated in Table [I] which shows the values of vari- 
ous contributions to the self-energy for a range of reduced 
temperatures. As in the homogeneous case, the removal 
of the infrared divergence can be seen in the substantial 
cancellation between the total static shift AE 1 -^ and the 
dynamic contribution AE^ D > (u>), while the small remain- 
ing difference AE (of order a tenth of a trap frequency 
or less) represents the overall change in energy compared 
to the Bogoliubov theory. 

Another interesting feature of the results in Table |U is 
the enormous size of the shifts AE4 and AE\ at high 
temperature. This is actually a separate issue from the 
question of infrared divergence [both these terms are part 
of the overall static shift AE^ S '] and has its origin in the 
large single-particle contribution to the non-condensate 
density n(r) which exists at finite temperature. This is 
proportional to (a positive power of) the reduced tem- 
perature t = T/T® rather than particle interactions and 
hence is large for t ~ 0.9. However, it is shown in Ref. Q 
that this single-particle contribution mostly cancels be- 
tween AE4 and AE\, as can be seen from their differing 
signs in Table U It is the residual large size of the sum 
A £4 + AE\ sa AE^ which represents the infrared di- 
vergence problem and which is removed by the inclusion 
of the dynamic term A£' D '(w). Further discussion of 



this point can be found in Ref. 0- 

Although the substantial cancellation between the 
static and dynamic terms means that convergence in the 
total self-energy is much better than in these individual 
contributions, it also means that small fractional errors 
in any one term can translate to large errors in the final 
answer, so care is required to preserve high numerical 
accuracy. A further difficulty in this regard comes from 
the fact that a number of quantities in the theory typi- 
cally converge slowly as the numerical cutoff energy E cut 
is increased. A simple example occurs in the calcula- 
tion of the total number of non-condensed atoms N nc , 
where at high temperature most of the particles reside in 
non-degenerate single-particle states above the numerical 
cutoff. If such quantities are required (as they may be; 
for example we need to calculate N nc in order to calcu- 
late the reduced temperature t) it is essential to include 
the effect of high-energy levels above the cutoff using the 
semi-classical approximation described in Appendix [BJ 

The semi-classical approximation can also be used to 
improve the convergence of our numerical results for the 
self-energy. Once E cut is greater than of order \Qfuo r 
this approximation is highly accurate for all equilibrium 
quantities, and hence also for the static shifts. Indeed if 
it is only required to calculate static terms (as in the var- 
ious versions of the HFB theories for example [3^ ) then 
rather small values of E cut can be used while maintain- 
ing high accuracy. Unfortunately, the semi-classical ap- 
proximation for the dynamic shift AEp D ^ (u>) is much less 
accurate (at least in the simplified form in which we have 
implemented it, cf. Appendix IB"|) . This is awkward be- 
cause we are ultimately interested in the small difference 
between the static and dynamic shifts and must there- 
fore calculate these terms to the same level of accuracy. 
One way to overcome this problem is by brute force, and 
we have therefore used a very large value for E cut in our 
simulations. This reduces the importance of the semi- 
classical terms while increasing their accuracy to roughly 
5 — 10% for the dynamic self-energy. Our final results 
are therefore converged to within about 10~ 2 frui r at the 
highest temperatures we consider. 

However, we have also developed a more sophisticated 
and vastly more accurate solution to this problem which 
takes into account the different ways in which a numerical 
cutoff should be introduced into the summations defining 
the static and dynamic terms. All static terms involve 
equilibrium non-condensate mean-fields and are therefore 
defined by a single summation over quasiparticle modes 
with a summation label i [cf. the time-independent limit 
of Eqs. ©-©I- Numerically, this summation is carried 
out for all quasiparticles modes with energies below the 
cutoff, ej < E cut . In contrast, the dynamic terms involve 
fluctuations in the non-condensate mean-fields and are 
calculated from a double summation over quasiparticle 
modes with summation labels i and j as in Eqs. (|A2|) - 
(|A4|) . These expressions are derived by allowing the 
quasiparticle wavefunctions Uj(r, t) and Vi(r, t) to fluc- 
tuate by writing Ui(r, t) = Ui(r) + Sui(r,t) etc and ex- 
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TABLE I: Contributions to the self-energy for the m = 2 mode and a range of reduced temperatures t for the parameters of the 
JILA TOP trap and a condensate population of Ao = 6000. All terms include a semi-classical contribution from high energy 
states, UV renormalization (if appropriate) and are given to two decimal places in units of huj r . The various contributions to 
the static shift are introduced in Eq. 1261 . the dynamic shift A£ (D) H is evaluated at the unperturbed Bogoliubov frequency 
tp/h, AEq is the contribution to AE^ D \uj) from the zero mode [cf. Eqs. 1A12I and HA13H and the final column AE is the 
total energy shift obtained as the sum of AE^ S ^ and AE^ D \ The results for the dynamic shift are calculated using a symmetric 
summation and are therefore too large by about 0.01 (see text and Fig. 0. The values of AE in the final column should 
therefore be reduced by this amount, with the result that the zero temperature shift is essentially zero |3S)l . 



panding the fluctuations using the static quasiparticles 
as a basis via Sui(r, t) — J2j Xij(t)uj(r) as in Eq. i|Al(l . 
This expansion is the origin of the second quasiparticle 
index j in the dynamic summations. 

The two summation indices therefore have different 
origins and should not be subject to the same numeri- 
cal cutoff. The index i must encompass the same states 
in both the static and dynamic terms as it ultimately 
describes which quasiparticle modes are included in the 
definition of the non-condensate mean-fields. The basis 
associated with the index j must be able to describe the 
dynamics of all these modes correctly, especially those 
with energies near E cut , and must therefore include a 
greater range of states since the dynamics of a mode near 
the cutoff may have a significant overlap with states of 
higher energy. We must therefore associate a different nu- 
merical cutoff with the two summation indices and take 

(?) (i) 

E^Jt > E c ^t- This asymmetric summation ensures that 
the static and dynamic terms are calculated with com- 
parable accuracy and produces a controlled cancellation 
between the corresponding self-energies for a finite ba- 
sis. For the results presented in this paper we have taken 

(i) (i) 

E^rJt = E c ^ t + 20hw r which is sufficient to provide the 
convergence we require. 

If we are only interested in the total self-energy (rather 
than individual contributions to it) the use of an asym- 
metric summation negates the need for the semi-classical 
approximation and leads to much more rapid and reli- 
able convergence, as shown in Fig. ^ We also see that 
this new technique leads to a small but non-negligible 
shift in the energies at all temperatures. In particular, 
the shift of order 0.01huj r seen at zero-temperature in 
the final column of Table [I] and in our earlier work |{| 
disappears if the dynamic self-energy is calculated using 
the new asymmetric summation method 38] . The new 
results are more consistent with analytical expectations 
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FIG. 1: Convergence of the quasiparticle energy E p defined 
in Eq. (13411 for the m — 2 mode as a function of numerical ba- 
sis energy cutoff for the parameters of the JILA experiment. 
Upper curves are at T — 0, lower curves are at T = 300nK 
(t ~ 0.88). Squares and solid line: asymmetric summation 
without the semi-classical approximation; circles and dashed 
line: symmetric summation including the semi-classical ap- 
proximation; asterisk with dot-dashed line: symmetric sum- 
mation without the semi-classical approximation. The lines 
are provided as guides for the eye. 



|17| and are more reliable. 

A difficulty with the new approach, however, is that 
the different ranges for the summations in the dynamic 
terms mean that in general the expressions for the sum- 
mands can not be symmetrized with respect to the labels 
i and j. This complicates the calculation, but it turns 
out that it is only an issue in practice for the zero tem- 
perature contributions (see Appendix |A"|) . Nonetheless, 
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for this reason, we have had to rederive the formulae for 
the dynamic terms being careful not to symmetrize with 
respect to the quasiparticle indices. The new expressions 
are given in Appendix ^ and replace the corresponding 
ones given in 3| where symmetrization was frequently 
used to simplify the calculation. Of course, the new re- 
sults reduce to those given earlier in the case that sym- 
metrization is allowed and in particular the expressions 
coincide in any exact calculation where both £^« t an< ^ 
E^ t are infinite. 

The substantial cancellation of the large static and dy- 
namic terms suggests that it may also be possible to avoid 
these issues by an appropriate reformulation of the the- 
ory. One possibility, first introduced by Popov, is to use 
density/phase variables for the perturbation theory with 
low-energy states. A similar approach has recently been 
developed and applied to quasi-condensates by Mora and 
Castin [3^| and it would be an interesting subject for fu- 
ture work to rephrase the current calculation using this 
formalism. 



F. Finite size effects 

The results in Table [I] provide interesting information 
on the importance of various finite size contributions to 
the energy shift. Interactions between the equilibrium 
non-condensate mean-fields and the condensate affect the 
static condensate shape and since this provides the mean- 
field for the quasiparticles there is a second order effect 
on their energy described by AE'shapo- This contribution 
is absent in the homogeneous limit and is therefore a 
finite size effect, but it is clear from the table that it can 
be very significant. Indeed for the particular parameters 
of the JILA experiment, this contribution is roughly 4-5 
times larger than the overall energy shifts. 

It is also interesting to see the importance of the con- 
tributions AE^ and AEq which arise from a careful 
treatment of wavefunction orthogonality in the number 
conserving approach. The function /(r,i) in the gener- 
alized GPE of Eq. Q arises specifically from the fact 
that the non-condensate is defined to be orthogonal to 
the condensate, while the contribution from AEq to the 
dynamic shift describes the explicit effect of the two zero- 
energy quasiparticle modes to the dynamics of the non- 
condensate [see Eqs. (jA12|) - (|A13f) ]. While both these con- 
tributions are small compared to the other shifts they 
are both comparable in size to the final answer and are 
therefore in principle significant. We should also point 
out, however, that AEj (p) does not give the full con- 
tribution of the function f(r,t) in the generalized GPE 
because there is also a dynamic contribution included in 
the quoted results for AE^^p). It is also clear from 
the table that there is substantial cancellation between 
AEj S \p) and AEo which mitigates their effect. How- 
ever, we have calculated the predictions of the SOBP 
theory neglecting the explicit effects of both f(r,t) and 



AEo (but continuing to use static quasiparticles orthog- 
onal to the condensate) and have found that the results 
for the energy of the in = 2 mode clearly disagree with 
the measurements obtained at JILA. This calculation and 
its implications will be reported and discussed elsewhere 

m 

This raises an interesting and important question for 
future work which is to what extent the effects of these 
terms are reproduced in a broken symmetry description 
in which the quasiparticle wavefunctions are not orthog- 
onal to the condensate, the function /(r, t) does not ap- 
pear and the quasiparticle description has to be supple- 
mented with the 'missing eigenvector' if it is to form a 
basis. We presume that a sufficiently careful treatment of 
these issues should be equivalent (for large condensates) 
to the results given here but the relative size of these 
terms in our calculation for the JILA experiment shows 
that such issues should not simply be ignored. 



G. Extracting energies and decay rates from the 
response function 

If the self-energy Y, p (l)-\- ij) and the driving of the con- 
densate by the thermal cloud [described by AP^^uj + 
ij) I P p q(oj)] are roughly independent of frequency, the 
energy shift can be calculated straightforwardly from the 
poles of G p (uj + ij) by finding the solutions to 

E p = tiuj p = Re [E p {uj p + , (34) 

with E p (ui) = e p + E p (a;) as in Eq. I|24|) . The correspond- 
ing decay rate is then given by 

Tp = -lxn[Ep{up + i 1 )]/h. (35) 

An example of this procedure is shown in Fig. [21 for the 
m = 2 mode in the JILA experiment. As can be seen 
the self-energy is relatively smooth near the unperturbed 
quasiparticle resonance but has significant structure fur- 
ther out. 

The situation of a smooth self-energy arises when an 
excitation couples to a continuum of decay channels, as 
in the homogeneous limit, and leads to resolvents Q p and 
IZp which are lorentzians. For a finite system, however, 
Y>p(uj + 17) depends on frequency, and both Q p and 1Z P 
can differ significantly from perfect lorentzians. In this 
case the line shape depends on the details of the system, 
and we have to extract energies and decay rates by fitting 
the response to a suitable function. When the spectrum 
is strongly distorted (as it can be at high temperature for 
example), the results obtained from these fits can be sen- 
sitive to the exact fitting procedure and this ultimately 
puts error bars on the theoretical predictions. We have 
therefore implemented two fitting procedures to deter- 
mine the energies and decay rates: in our earlier work 
we used a complex lorentzian to fit the full response in 
the frequency domain |5j, while in the present paper we 
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FIG. 2: Real and imaginary parts of the frequency-dependent 
quasiparticle energy E p (u) of Eq. 12 H for the m — 2 mode 
and reduced temperatures t = (top line), 0.65 (middle line), 
and 0.9 (bottom line), a) Im[E p (u>)], b) Re[E p (ijj)], both in 
units of huj r . In b) the diagonal is the line E = fuo and the 
intersections with the curves give the poles of the resolvent 



The corresponding values of T p can then 



E p of Eq. 

be read off from a) using Eq. 13511 . The vertical dashed line 
gives the position of the unperturbed quasiparticle resonance 
e p /huJr = 1.44 to 3SF. 



mimic the experimental procedure exactly and fit a de- 
caying sinusoid to the response in the time-domain. 

We find that at low temperatures both methods pro- 
duce excellent fits and the energies and decay rates can 
be extracted straightforwardly. At higher temperatures, 
however, the spectra develop noticeable non-lorentzian 
structure (see Figs. [3] and and it becomes more diffi- 
cult to obtain reliable fits. This is not surprising when 
one considers that a fit is an attempt to model a complex 
spectrum with very few parameters which is difficult to 
do reliably when the spectrum does not have the assumed 
form. The problem is exacerbated from the fact that the 
intrinsic width of the spectra at high temperature is of 
order Q.2huj r (cf. Fig.0) while we would like to determine 
the central frequency to an accuracy roughly an order of 
magnitude better than this. 

These problems can be overcome to some extent by 
including the spectrum of the perturbation P p q(uj) as 
a known weight function in the fit. This ensures that 
only the experimentally relevant range of frequencies are 
included and has the effect of suppressing some of the 
non-lorentzian structure in the wings of the distribution. 
Specifically, we fit the spectra in the frequency domain 
to the function 



(36) 



where the fit parameters are the real constants E„ and T p 
and the complex constants A and C [C accounts for any 
linear frequency dependence of AP^ (lu)]. The experi- 



mental resolution 7 is subtracted from the fitted width 
parameter T p to give the intrinsic width. The weight 
function P p q(uj) is calculated using the functional form 
given in Eq. (|33|) and has the frequency dependence 



P p o(lo) cx 



sm[(w - uj d )t d ] 

{uj - ljJ d )td 



(37) 



The results are not particularly sensitive to the choice of 
the central drive frequency u>d, which was chosen in the 
experiment to maximise the observed response. In our 
calculations we chose it to be equal to the unperturbed 
Bogoliubov energy e p at T = 0, while at higher temper- 
atures we take it to have the value of E p found at the 
closest lower temperature for which we have data. In 
this way the trend of the central drive frequency broadly 
follows the trend of the energy shift. 

For the time-domain fits, we take the Fourier transform 
of the response functions (including the factor P p o(uj) as 
above) and fit the real part (which describes the real con- 
densate density oscillations) to a decaying sinusoid of the 
form Ae~ Tpt sm(E p t+<j)) + B where the fit parameters are 
all real. The times t are chosen in the range td < i < t t, s 
where t Q b s is the experimental observation time, which is 
34ms in the JILA experiment Q . The resolution param- 
eter 7 (which is no longer required) is removed by multi- 
plying the oscillations by e +7 * before fitting. The results 
from this procedure are in complete agreement with the 
fits in the frequency domain at low temperatures where 
the spectra are well described by lorentzians, but at high 
temperatures there are differences of the order of a few 
\Q~ 2 hw r . This uncertainty in how to extract energies and 
decay rates from non-lorentzian spectra therefore repre- 
sents the theoretical uncertainty in the predictions. 

Looking ahead to the results shown in Fig. El we see 
that the energy E p extracted from assuming a frequency- 
independent self-energy is clearly in better agreement 
with experiment for the m — 2 mode at the highest tem- 
perature than the values extracted from the fits. Inspec- 
tion of the shape of the frequency-dependent quasiparti- 
cle energy for this case, given in Fig. [3 shows that the val- 
ues of E p are indeed representative of the typical values in 
the frequency range of interest. Given the arduous nature 
of the calculation, it is therefore somewhat frustrating 
that the agreement with experiment should be substan- 
tially worsened by the fitting procedure. Nonetheless, 
fitting the condensate density oscillations in the time- 
domain mimics the experimental procedure and is there- 
fore presumably the appropriate method to use when the 
spectrum has non-lorentzian structure. 



IV. NUMERICAL RESULTS 

In this section we present a detailed analysis of the 
SOBP predictions for the 1997 JILA experiment 3- In 
particular, we show the functional form of the response 
functions including and excluding the effect of direct ex- 
citation of the non-condensate and discuss the under- 
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lying physics. We then present results for the ener- 
gies and decay rates of the lowest energy to — 2 and 
to = modes, corresponding to the quantum numbers 
(m,p, n) = (2, 0, 1) and (0, 0, 1) respectively. In Sec. lTvTl 
we provide a detailed analysis of the experimental results 
for condensate population Nq and temperature T to de- 
termine the appropriate values to use in our simulations 
and the uncertainties which variations in these quanti- 
ties can be expected to produce in our results (with the 
exception of the analysis in this section we have consis- 
tently taken N = 6000). We also present results for 
the axial and radial dipole modes (quantum numbers 
(m,p, n) — (0, 1, 1) and (±1, 0, 1) respectively) and show 
that these are obtained correctly in the SOBP theory pro- 
vided that excitation of the non-condensate is included. 
Finally, we discuss the relative phase of the condensate 
and non-condensate density oscillations. 



A. Response functions 

In this section we consider the functional form of the 
resolvents G p {oj) and 1Z p (oj) defined in Eqs. and ij2*T|) 
respectively for the lowest energy to = 2 and m = 
modes. For the m = 2 mode, symmetry considera- 
tions mean that there are no fluctuations in the con- 
densate population (SNq(uj)—0) and these quantities are 
directly proportional to the condensate density response 
in the frequency domain. For the m = mode, how- 
ever, SNq(u>) is non-zero and should be included in the 
analysis of the density response. This can be done as 
described in Appendix IA 21 but has very little effect on 
the detailed results. For the purposes of the present dis- 
cussion, the difference is not substantial enough to merit 
complicating the issue further and we will therefore treat 
G p (uj) and TZ p (uj) as giving the density response directly 
for both modes. The only exception to this comes in 
Sec. IIV El where we consider the relative phase of the 
condensate-non-condensate oscillations and include the 
effect of &Nq(<jj) for the m = mode. 

The absolute values of these response functions are 
plotted in Figs. |3| and 0] for the to = 2 and to = modes 
respectively as functions of frequency uj for a range of 
reduced temperatures. In each case, the upper panel (a) 
shows Q p (u)), which neglects direct driving of the non- 
condensate by the perturbation, while the middle panel 
(b) shows the full response 1Z p (uj) where this process is 
included. The lower panel (c) shows the contribution to 
the self-energy from a subset of low-energy Landau and 
Beliaev processes and is included to aid analysis of the 
response function plots. 

We focus first on the results for Q p (lj). As can be seen, 
at t = the response of both the to = 2 and to = 
modes is a lorentzian positioned almost exactly at the 
frequency of the unperturbed Bogoliubov energy In fact, 
we find that the shift in the frequencies at zero temper- 
ature is essentially zero for both modes (see Fig. [JJ, in 
contrast with our earlier calculations where shifts of or- 




FIG. 3: a) 7 x |£/j>(w)|, b) 7 x \1Z v (lo)\: modulus of resolvent 
as a function of frequency for the m = 2 mode [(m,p,n) = 
(2,0,1)] and t = (solid), t = 0.65 (dotted, x2), and t = 
0.9 (dashed, x3). c) — Im[£ p (w)] at t = 0.9 (dashed curve) 
and a few contributing Landau and Beliaev processes (solid 
and thick-dashed vertical lines respectively) . These processes 
are drawn at their resonance frequencies Wy with a height 
corresponding to their amplitude in the self-energy. The three 
largest contributions are labelled with the (m,p,n) quantum 
numbers for the two quasiparticle modes involved (the top 
label applies to the mode with the lower energy). The thin, 
vertical, dashed line indicates the position of the unperturbed 
Bogoliubov resonance, e p = 1.44fiw r to 3SF. 



der 2 x I0~ 2 huj r were found 51. The new results are 
consistent with the analytical results of Giorgini 01 ob- 
tained in the Thomas-Fermi regime, and are expected 
to be more accurate because they incorporate the asym- 
metric summation method described in Sec. IIII El which 
has a substantial effect at zero-temperature, as shown in 
Figffl The width of the response functions at zero tem- 
perature is entirely due to the parameter 7 introduced 
to model the finite experimental resolution. As the tem- 
perature increases, however, the response develops an in- 
trinsic width and there is a noticeable downwards shift 
in the position of the peak with a concomitant decrease 
in height. 

At the highest temperatures, considerable structure is 
observable in the wings of the response, which is no longer 
a simple lorentzian. The physical processes responsible 
for this structure can be identified using the self-energy 



13 




FIG. 4: Plots as in Fig.^Jbut for the m = mode [(m,p, n) — 
(0,0, 1)]. The curves in a) and b) now correspond to: t — 
(solid), t = 0.65 (dotted, xl.5), and t = 0.9 (dashed, x2). 
The unperturbed Bogoliubov resonance is at e p — 1.85/iuv to 
3SF. 



plots of Figs. |3t and Here we show the negative of 
the imaginary part of the self-energy at a high temper- 
ature t = 0.9 (the dashed curve) as well as vertical bars 
whose height and position indicate the contribution from 
individual Landau and Beliaev processes. As Eq. I|27(l 
demonstrates, the imaginary part of a dynamic term is 
formed from the combined effect of many such processes, 
each broadened by the resolution parameter 7. Where a 
large number of these resonances overlap, the self-energy 
is smooth but if a few large resonances dominate at some 
frequency then sharp features are observed leading to 
non-lorentzian structure in the resolvent. This structure 
generally occurs at a frequency slightly further from the 
Bogoliubov energy e p than the features in the self-energy 
because of level repulsion (typical of a second order per- 
turbation calculation) between these processes and the 
central resonance. 

For the to = 2 mode, Figs. |2t and |3J; show that the 
structure in G p (u>) around u /u) r = 0.7 is a consequence of 
a few Landau processes of the form (2, 0, 1) + (to, 0, n) — > 
(to + 2,0, n) for small to, while that near u>/u> r ~ 2 
is a result of a single strong Beliaev decay into two 
x-y dipole modes, (2,0,1) -> (1,0,1) + (1,0,1). For 
the to = mode, the structure around oj/uj t ~ 1.3 in 
Fig. is the result of a few Landau processes of the 



form (0, 0,1) + (m, p, n) — > (to, p, n+ 1) for small values of 
to. For both modes, the processes described involve low- 
energy collective excitations of the thermal cloud, with 
the result that they have quite large matrix elements and 
hence can be significant at finite temperature when the 
associated rates are strongly enhanced by Bose stimula- 
tion. For the TOP trap geometry, these processes oc- 
cur at frequencies quite a long way from the principal 
resonance, but the large width of the response at finite 
temperature means that they can still cause a significant 
distortion of the spectrum, even to the extent that the 
greatest response no longer occurs in the vicinity of the 
original Bogoliubov mode. The effect of these processes 
could be increased by tuning the trap geometry to shift 
them closer to resonance, although a full discussion of 
their importance requires the inclusion of direct excita- 
tion of the non-condensate by the perturbation, which 
we now consider. 

The effect of including direct driving of the non- 
condensate is shown in Figs. |3Jd and QJa, where the full 
response function 7Z p (u>) is plotted. For the to = 2 mode 
this process only affects the wings of the response, which 
is otherwise qualitatively the same as if it is neglected. 
The main change is an enhancement of the structure 
around u> /u) r ~ 2, which is again a consequence of a single 
Beliaev process in the thermal cloud involving the exci- 
tation of two (1, 0, 1) dipole modes. For the m — mode, 
however, there is a dramatic change in the form of the 
response function at high temperature, with a growing 
peak at uj/u) r = 2 which eventually dominates the spec- 
trum. In this case, the change in the response is due, not 
to a single Beliaev process, but rather to a large number 
of weak Landau processes. The external perturbation is 
proportional to r 2 and hence couples strongly to high- 
energy, single-particle modes in the thermal cloud. Since 
these modes are weakly-interacting, there are a large 
number of Landau processes with frequency differences 
near the ideal gas value of 2uj r . When the temperature 
is high enough that these modes are significantly pop- 
ulated, the non-condensate has a large response at this 
frequency, which it can transfer to the condensate via 
their dynamical coupling. This process ultimately dom- 
inates the condensate response, with the result that at 
high temperature it is more likely to be excited indirectly 
via the thermal cloud than directly via the perturbation. 
This provides the microscopic explanation for the strong 
upward shift in the excitation energy of the to = mode 
at t ~ 0.6 observed in the 1997 JILA experiment (see 
Fig.© 0,0. 

The importance of including the direct excitation 
of the non-condensate by the external perturbation is 
perhaps shown most clearly in Fig. [S] where we plot 
Qp{oS)/Q p {oS) which describes the ratio of the indirect ex- 
citation of the condensate via the thermal cloud to its di- 
rect excitation by the external perturbation [cf. Eq. (|22() ]. 
At low temperatures and for most frequencies, this ra- 
tio is much less than one indicating that the condensate 
is mostly driven by the external perturbation. However, 
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FIG. 5: Plot of the absolute value of Q p {ui)/Q p (uj) as a func- 
tion of frequency for the m = mode. As shown by Eq. 1221 . 
this gives the ratio of indirect condensate excitation via the 
thermal cloud to its direct excitation by the external pertur- 
bation. Curves are shown for t = (solid), t — 0.65 (dotted), 
and t = 0.9 (dashed) as in Fig. 0] 



there is also a peak at to = 2u r which grows with temper- 
ature and at t = 0.9 we find that the condensate is driven 
roughly five times as strongly via the thermal cloud as it 
is directly via the perturbation. 

We should note that whereas the results for G p (cu) are 
fundamental, depending on the intrinsic couplings be- 
tween the condensate and thermal cloud, the response 
function TZ p (u>) depends on the assumed form of the ex- 
ternal perturbation and as such will differ in detail from 
one experiment to another. This provides a handle to ex- 
plore the relative importance of the two mechanisms for 
exciting the condensate. For example, if the perturba- 
tion is chosen to be spatially localized around the centre 
of the trap (rather than the r 2 form assumed here) then 
the effect of excitation of the thermal cloud will be greatly 
reduced and the condensate response should be well de- 
scribed by Qpiui) alone. It would be interesting therefore 
if the response of the m = mode was remeasured us- 
ing perturbations of different spatial form so that the 
downward shift for this mode predicted in the absence of 
thermal cloud driving can also be observed (see Fig.EJ). 



FIG. 6: Ab initio theoretical excitation energies E (open sym- 
bols) compared with experiment (filled circles with error bars) 
for (a) the m = and (b) the m = 2 mode. Diamonds neglect 
direct thermal driving (Q p ), open circles include it (1Z P ) and 
squares give E p of Eq. (13411 . The dashed line is the Bogoli- 
ubov energy e p . Differences between diamonds and squares 
are due to non-Lorentzian structure in Q p . There are no free 
parameters in the theoretical results. 
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FIG. 7: Theoretical decay rates (r) compared with exper- 
iment for (a) the m = mode and (b) the m = 2 mode. 
Symbols as in Fig. 



B. Energy shifts and decay rates 

We extract energies and decay rates from the resolvents 
given in Figs. and 21 by finding poles of the self-energy 
and by using fits to the oscillations in the time-domain, as 
described in Sec. IIII Gl The results of these calculations 
for the m — and m = 2 modes are shown in Figs. Eland 
[7| and are essentially the same as those given in Ref . [5| . 
There are two minor differences compared with our ear- 
lier work, however. The first comes from the use of the 



asymmetric summation method which shifts the energies 
down slightly at all temperatures and reduces the zero- 
temperature shifts to almost zero. The second comes 
from fitting the spectra in the time-domain rather than 
the frequency domain which makes the upward shift of 
the m = mode at t ~ 0.6 somewhat less abrupt when 
direct driving of the thermal cloud is included. However, 
all our earlier conclusions remain unchanged, and as we 
commented in Sec. lIII (Ti the difficulty in extracting mean- 
ingful fits from non-lorentzian spectra represents the the- 
oretical uncertainty in the final predictions. 
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It is perhaps worth noting that a critical analysis of 
Fig. |nt> suggests that the theory slightly underestimates 
the shifts seen in the experiment at high temperature for 
the m = 2 mode. Although the error bar on the final 
point is rather large in this case, this may be an indi- 
cation that the experiment is seeing effects beyond the 
SOBP theory. Recently, Liu et al. 40] have included 
self-consistent effects of the non-condensate mean-field 
dynamics for a spherical trap geometry and shown that 
this has the effect of enhancing the shifts predicted by the 
SOBP theory for the monopole mode. It would therefore 
be interesting to apply this approach to the JILA exper- 
iment to see if it improves agreement with experiment 
for the m = 2 mode at high temperature. However, 
the calculation of [4GJ neglects collisional effects in the 
non-condensate which may also be significant in a sys- 
tematic extension of the SOBP theory to next order in 
the dilute gas parameter. These effects are included in 
the recent theories of Walser et al. |4lL |42| , Wachter et 
al. 0] and Proukakis 0], although to the best of our 
knowledge these approaches have also not been applied 
to anisotropic geometries. 



C. Effect of changes in condensate population 

For a fixed trap geometry and atomic species, the main 
input parameters to the numerical calculations are the 
condensate population No and the absolute temperature 
T in nanokelvin. The relevant experimental data for 
these quantities is given in Fig. 1 of Ref. Q and is re- 
produced in Fig. [SJ where Nq, T and the total atom 
number TV are plotted as functions of reduced temper- 
ature t 0]. In order to verify the correct inputs for our 
calculation, we have used the ordinary Bogoliubov the- 
ory with Nq and T as input to calculate the total atom 
number N from the number of non-condensed atoms 
N nc (T, N ) = Jdrn(r) via N(T, N ) =N + N nc (T, N ). 
This is then used to calculate the reduced temperature t 
using Eq. (|2*9")l to obtain the ideal gas critical tempera- 
ture T®. The results of these calculations for four values 
of No and a range of temperatures are shown in Fig. [S] as 
solid lines. 

This calculation shows that while the experimental 
data for Ao is consistent with a condensate population 
in the region of Ao = 6000 or greater for t < 0.9, the 
results for N and T are more consistent with No — 3000. 
Since the Bogoliubov theory should adequately describe 
the thermodynamics of the experiment, we conclude that 
there may be substantial error (possibly systematic) in 
some of the experimental data. We stress that these nu- 
merical calculations come from a full implementation of 
the ordinary Bogoliubov theory for the anisotropic trap 
geometry and finite particle number of the experiment 
and the convergence on all numerical quantities shown is 
of order a part in 10 5 . 

We have assumed that the data for the condensate 
population is probably the most reliable and for this 
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FIG. 8: Bogoliubov theory (solid lines) compared with exper- 
iment (open circles) for (a) T in nanokelvin vs t, (b) N vs t, 
and (b) N vs t. Solid lines correspond to N = 8000, 6000, 
3000, and 1000 from top to bottom in each case. 



reason we have used a value of Ao = 6000 for most of 
our calculations and for the results reported in Ref. |5j]. 
Given the uncertainty in Ao, however, it is important 
to calculate what effect this has on our earlier predic- 
tions. In Fig. El we plot the energies obtained from the 
SOBP theory for both the m = 2 and m — modes for 
N = 1000, 3000, 6000 and 8000, which cover the range 
of possible values in the experiment. As can be seen 
the results are relatively insensitive to the value of A , 
although the case Ao = 1000 is clearly excluded. In- 
deed, at high temperature the difficulty in extracting a 
meaningful energy from a non-lorentzian spectrum leads 
to a similar uncertainty in the prediction as variation of 
the condensate population within the relevant range. We 
conclude that the uncertainty in the relevant values of Ao 
and T to use in our simulations has little effect on our 
overall results, although clearly it would be advantageous 
to have new experimental results where this uncertainty 
was removed. 
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FIG. 9: Quasiparticle energies for a) the m = mode and 
b) the m — 2 mode as a function of reduced temperature t 
for No = 1000 (triangles), iV = 3000 (stars), 6000 (squares) 
and 8000 (circles). The results for m = and the upper set 
of points for m — 2 come from fits to the resolvent 1Z P in 
the time domain, while the lower set of results for m = 2 
correspond to the value of E p from Eq. (1341 . The difficulty in 
extracting meaningful energies from non-lorentzian spectra is 
visible by the large difference between these two sets of results 
for t > 0.8. 



D. Dipole modes 

In a trapped system at finite temperature both the con- 
densed and non-condensed atoms can undergo centre-of- 
mass oscillations, corresponding to the excitation of the 
dipole modes of the system. If there is relative motion of 
the two clouds then these oscillations will ultimately be 
damped 6J. In a harmonic potential, however, the gener- 
alized Kohn theorem [2^| shows that in-phase oscillations 
are an exact solution of the full equations of motion and 
are completely decoupled from any internal dynamics of 
the clouds, independent of both temperature and par- 
ticle interactions. As a consequence we have the exact 
result that the system has undamped dipole oscillations 
with energies corresponding exactly to the principal trap 
frequencies. Whether or not these modes are obtained 
correctly is therefore an important test of any theoretical 
description of a condensed system at finite temperature. 

In Ref. we showed that the SOBP theory is con- 
sistent with the generalized Kohn theorem for the dipole 
modes to within the small parameter of the theory. This 
is only the case, however, if the effect of the external per- 
turbation on the non-condensate dynamics is included. 
This makes sense because it is clearly not possible to 
describe an in-phase oscillation of the condensate and 
non-condensate if the force which generates the motion 
is assumed at the outset to act only on the condensate. 

Our numerical results for the axial dipole oscillation 
are shown in Fig. llUI where we plot the magnitude of the 
response functions Q p {uS), Q p {lo) and 7Z p (u>) for a range of 
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FIG. 10: Plot of the absolute value of the resolvents (a) Q p (lu), 
(b) G v (u->) and (c) lZ p (u>) as a function of angular frequency ui 
for the dipole oscillation along the 2-axis, (m,p, n) — (0, 1, 1). 
In all cases the spectra are shown for reduced temperatures 
t = 0,0.15,0.3,0.45,0.6,0.75, and 0.9. In (a), temperature 
increases from top to bottom, in (b) it increases from bottom 
to top, while in (c) the various curves are indistinguishable 
on this scale. 



temperatures. In Fig. HOfa . the effect of direct excitation 
of the non-condensate is excluded and so the dipole mode 
is not obtained exactly. In this case, the perturbation 
moves the condensate relative to the centre of the non- 
condensed atoms with the result that the oscillations are 
damped. This is visible as an increase in the width of the 
response function and a decrease in its peak height as the 
temperature increases (there is no significant change in 
the position of the peak frequency) . The damping is very 
light compared to the other low-energy modes, however: 
at t ~ 0.9 we have T p ~ 0.03w r , which is 20% or less of 
the decay rate for the m — and m = 2 modes shown in 
Fig.0 

It is also interesting to consider the condensate re- 
sponse in the case where it is driven only via the ther- 
mal cloud, corresponding to the experiment of Ref. [||. 
This is described by the response function Q p (u>) defined 
in Eq. ill'l'li and is shown in Fig. 110b for the case that 
P(r, ui) oc z, which is the linear form required to excite a 
pure axial dipole oscillation. In this case the amplitude of 
the response increases with temperature as the size and 
importance of the thermal cloud increases. Although not 
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obvious from the figure, the width of this response func- 
tion is also slightly narrower than the parameter 7 we 
introduced into the resolvents to model the experimental 
resolution. This means that when we transform Q p (uS) 
to the time-domain and remove this (artificial) decay by 
multiplying by e +7 * we see a clear growth in the am- 
plitude of the condensate density oscillations with time. 
This is in contrast to the case for Fig. 110b , where the os- 
cillations decay with time and conforms precisely with 
our expectations: In this case the perturbation causes 
the thermal cloud to oscillate through the condensate re- 
sulting in a transfer of energy and momentum and hence 
to condensate oscillations which initially grow with time. 



The full condensate density response is obtained by 
adding Q p {oS) and Q p (oS) to obtain lZ p (u>), which is shown 
in Fig. II Oh . In this case the curves for the various tem- 
peratures collapse onto each other, and to high accuracy 
the response is a pure lorentzian centred on the axial 
trap frequency. As expected, there is no damping of the 
oscillation and the width of the response shown simply 
corresponds to our resolution parameter 7. Extracting 
the energy and decay rate from the spectra, we find that 
the peak of the lorentzian is at the axial trap frequency 
lo z to within 0.1% for all temperatures, which is of the 
order of our numerical accuracy, while the intrinsic decay 
rate is indistinguishable from zero. This result confirms 
that this dipole mode is obtained correctly in the SOBP 
theory and acts as an important check on our numerical 
method. 



We have also calculated the spectra for the dipole os- 
cillations in the x-y plane [quantum numbers (m,p, n) — 
(±1,0, 1)] with very similar results. If direct excitation 
of the non-condensate is included, we find that the mode 
frequency is obtained correctly to within 0.1% for t < 0.6, 
although the error rises to nearly 1% at t = 0.9. The in- 
creased error in this case may be due to the approxima- 
tion that only the positive frequency pole of the conden- 
sate response is relevant in the calculation 0\ . In reality, 
the full response also has a contribution from a pole at 
u) = —lo t and this is a factor of y8 closer to the positive 
frequencies of interest for radial oscillations than it is for 
axial oscillations in a TOP trap geometry. 



Finally, we note that the use of the asymmetric sum- 
mation technique described in Sec . IIII El greatly improves 
the accuracy of our calculation of the dipole modes, es- 
pecially at zero temperature. If a symmetric summation 
is used instead, we find a systematic upwards shift in the 
zero-temperature frequency of both the radial and ax- 
ial dipole modes of order 2 x 10~ 2 clV, representing a 2% 
error for the radial modes. The asymmetric summation 
reduces this error by an order of magnitude to within the 
level of our numerical accuracy. 



E. Relative phase of condensate—non-condensate 
oscillations 

An important issue in the study of condensate dynam- 
ics at finite temperature is the relative phase of the oscil- 
lations of the condensed and non-condensed atoms. In- 
deed, it has been argued byBijlsma and Stoof an d 
by Al Khawaja and Stoof [19( that the JILA results for 
the m = mode can be explained by a transition from an 
out-of-phase motion of the two clouds to an in-phase mo- 
tion at high temperature. We have found that the SOBP 
calculation essentially confirms this conclusion. Combin- 
ing the expression for the condensate density fluctuations 
in Eq. Q16[) (including the explicit effect of the condensate 
number fluctuations SNq(uj) for the m = mode) with 
that for the non-condensate in Eq. 1|A2(1 we can straight- 
forwardly calculate their relative phase. Of course, both 
these quantities are functions of position and frequency 
as well as temperature, so we first take a suitable moment 
of both oscillations as described in Eqs. (|A20|) and (|A21|) 
and then evaluate their relative phase at the frequency 
which gives the best fit to the condensate oscillations (av- 
eraged over a region of 27) . This gives a measure of the 
relative phase of the two oscillations at the peak conden- 
sate response as a function of temperature. 

The results of this calculation for the m = 2 and m = 
modes as well as the axial dipole mode are shown in 
Fig. ^2 As can be seen, in all cases the oscillations are in 
phase near t = 0. If direct driving of the non-condensate 
is neglected, the oscillations become increasingly out-of- 
phase as the temperature rises, especially for the m = 2 
and m = modes. The results are very different when 
the effect of the perturbation on the non-condensate is 
included. In this case we see that the m = mode shows 
a distinct crossover to an in-phase motion at high tem- 
peratures. This is indeed what we would expect from our 
earlier analysis as in this case the condensate is driven 
predominately by the thermal cloud at a frequency above 
its natural Bogoliubov energy resulting in an in-phase re- 
sponse. A similar effect is also seen for the m = 2 mode 
although it is less strong because the thermal cloud res- 
onance at to = 2uj r is much less significant in this case. 
It is also gratifying to see that the condensate and non- 
condensate oscillations for the dipole mode remain locked 
in-phase at all temperatures when direct driving of the 
thermal cloud is included. This is consistent with the 
Kohn theorem and with the results we obtained for the 
dipole modes in Sec. lIVDl 



V. CONCLUSIONS 

In conclusion, we have presented a detailed analysis 
of the measurements of condensate excitations made at 
JILA in 1997 using a theory of the linear response of 
Bose-Einstein condensates at finite temperature that we 
have recently derived. We have shown the importance 
of including the direct effect of the external perturbation 
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FIG. 11: Plot of the phase of the non-condensate density os- 
cillations relative to those of the condensate as a function of 
reduced temperature for the parameters of the JILA experi- 
ment. The solid line gives the results for the m = mode, 
the dashed line for the m = 2 mode and the dashed-dotted 
line for the axial dipole mode. In each case the upper curve 
neglects direct driving of the non-condensate by the external 
perturbation while this effect is included in the lower curve. 



on the non-condensate dynamics and in particular we 
have demonstrated unambiguously that this provides the 
explanation for the anomalous behaviour of the m = 
mode observed at JILA and that it is necessary for a 
correct description of the dipole modes. A major issue 
in the numerical implementation of the theory is achiev- 
ing convergence of the results and we describe a novel 
asymmetric summation scheme we have developed which 
solves this problem. This makes it feasible to apply our 
theoretical formalism to a wide variety of more recent 
experiments involving larger condensate populations and 
more highly anisotropic traps. 



Sec. IIII El as a means of substantially improving the nu- 
merical convergence. The results given here replace the 
equivalent ones obtained in Ref. Q, to which they re- 
duce in the case of an exact calculation or a symmetric 
summation. We will, however, make frequent reference 
to this earlier work for the definition of numerous ma- 
trix elements which appear unchanged in the formulae 
derived here. 

Expressions for the fluctuations of all non-condensate 
mean-fields can be obtained by writing them in terms of 
time-dependent quasiparticle wavefunctions Ui(r,t) and 
Vi(r,t) and linearising these around their static values, 
Ui{v,t) — > Ui(r) + Sui(r,t) etc. These quasiparticle fluc- 
tuations are found by introducing an expansion in the 
equilibrium basis of the form 

(M-p^W < ai » 

and obtaining the expansion coefficients Xij from the 
solution of a linearized form of the time-dependent BdG 
equations of Eq. Q in the frequency domain. The results 
are given in Eq. (A33) of Ref. Q (see also Eqs. (IXfl-llXflll 
below) and can be taken over unchanged to the present 
work (the formulae for to not depend on whether we 
subsequently use a symmetric or asymmetric summation 
over the indices i and j). 

Following this procedure, the fluctuations in the non- 
condensate density and anomalous average are given by 

5n(r,r» = ^ [N lV *u* + (N + \)v*u*]X^(u) 

ij>0 
ij>Q 

+ WiUjU* + (Ni + lK*u,]A^H 

ij>0 

- Y WiUiU* + (Ni + ltyvAXjitu), (A2) 

ij>0 
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APPENDIX A: DEFINITION OF DYNAMIC 
TERMS WITH ASYMMETRIC SUMMATIONS 

In this Appendix we provide expressions for all dy- 
namic terms in the theory. These expressions are ob- 
tained without symmetrizing with respect to the quasi- 
particle labels involved and are therefore appropriate for 
the new asymmetric summation technique introduced in 



5m R (r,r» = ^ [N z v*v* + (N + 

ij>0 

+ Y WiUiUj + (N + l)u 3 u l ]X*_ :j (-uj) 

ij>0 

+ ]T [N iUj v* + (N + 1)v* U] ]X 13 (uj) 

ij>0 

- WiUiVj + (Ni + l)v*Ui\Xji(uj) 

ij>0 

+ N ^-$ 6$(r,uj)S(r-r'), (A3) 

where we have used the convention that the first wave- 
function has the spatial argument r and the second r'. 
The last line of Eq. (|A3|) is the contribution from the UV- 
renormalization. The summations in these equations are 
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over all positive-norm modes and therefore include the 
corresponding zero-energy mode with a population fac- 
tor N Q+ = (N.B. N 0+ ^ N ). 

All dynamic terms required in the theory are pro- 
jections of linear combinations of 5n(r,r',u>) and 
(5m R (r, r', u>), usually, but not always, with r = r'. A 
generic dynamic term can therefore be written in the 
form 



y>0 



2^ 
(A4) 



where the C P ij coefficients are matrix elements involv- 
ing integrals of quasiparticie wavefunctions whose defini- 
tion depends on the particular combination of 8h(r, r', u>) 
and 5m R (r, r', u>) appearing in the dynamic term in ques- 
tion. The population factor is equal to 2Ni for 
finite temperature contributions and unity for zero tem- 
perature contributions, while the final term A R (u>) is the 
UV-renormalization (if any). 

In an exact calculation where Y]^ ranges over all pairs 
of states, we can symmetrize the summands in the above 
equations with respect to i and j. Eqs. (| A2|) and (|A3|) 
then reduce to Eqs. (A26) and (A27) given in 0, and 
all expressions derived from them similarly reduce to the 
corresponding expressions given in Ref. [3j. In this case 
the Beliaev and Landau terms [respectively the first and 
second lines of Eq. <|A4fl ] acquire the familiar population 
factors of 1 + iV, + Nj and Nj — Ni respectively. However, 
any numerical calculation necessarily involves a finite ba- 
sis and, as discussed in the main text, we should take the 
associated cutoff at a higher energy for the j summation 
than for the i summation to improve the convergence 
properties. Thus numerically, the summation is 



(A5) 



ij>0 ij>0 



with E% > Ei and 0(x) the unit step function. This 
asymmetric summation requires the use of the expres- 
sions given in this appendix in place of those in Ref. • It 
also has the unfortunate side-effect that for many quanti- 
ties of interest, the matrix elements at zero-temperature 
differ slightly from their finite temperature values. This 
difference arises from the fact that the quasiparticie ex- 
pansion of h(r, t) has a different form at zero and finite 
temperature, cf. Eq. ©. 

The coefficients Xij in Eq. (|A1|) contain two distinct 
contributions, one describing excitation of the thermal 
cloud by the condensate and the other its direct exci- 
tation by the external perturbation, as in Eq. (A33) of 

, which we denote here by X^' and X-p respectively. 

(c) 

In general, X^ depends on a sum over all excited con- 
densate excitations, but if a single mode dominates the 



summation we can write X^ as 

X 13 (lo) = x[f{p,u)b p {u) + X\p{u), (A6) 

where b p (uj) is the condensate expansion coefficient of 
Eq. I|17|) • For the case that neither of the subscripts i or j 
refers to the positive-norm, zero-energy mode (i,j =/= 0+), 
Eq. (A33) of Ref. Q gives for the condensate contribution 



Y 



(B2) 



pi.) 



m~ 



-Y 



e, + ej 

(A) 



huj - 



Y v 

vn 



ei + e 3 
(si) 



Hu) — e, — 



while the contribution from the perturbation is 



x 



x 



(P) 



(CO) 



hu — a + ej ' 



hu 



X^J(-u) = 



hu — 



(A7a) 
(A7b) 
(A7c) 

(A8a) 
(A8b) 
(A8c) 



The matrix elements Y pi j, P^(u) and P£(uj) in these 
equations are defined in Eqs. (106), (113) and (114) of 
|3j. Their calculation involves integrals of three quasi- 
particie wavefunctions and the condensate (the Y pi j) or 
two quasiparticie wavefunctions and the external pertur- 
bation [P9(w) and 

If, either of the subscripts i or j refers to the zero- 
mode, however, the formulae are different. In this case 



all the coefficients x\p(u) are zero and the X$'(p,u>) 
are given by (cf. Eqs. (A35)-(A37) of 0) 



-X {c) *(n 



*Si(P,<") 

x[f_{p,-u) = -V ip , 



-u) 



i(P»w) = -w ip , 

(A9a) 

= U ip , (A9b) 



— n 



(A9c) 



where Uij, Vij and Wij involve integrals of products 
of two quasiparticie wavefunctions and are defined in 
Eqs. (99)-(101) of 0. 

We deal with the different expressions for these two 
cases by restricting the summation in Eq. (|A4|) to positive 
energy modes and calculating the zero-mode contribution 
separately. Our final expression for a generic dynamic 
term is therefore 



(A10) 



where the condensate contribution A p c J (u) is further sub- 
divided as 



aWm = a^h + a (p) + a r (p). 



pp 



(All) 
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Here A p c p +) (uj) and A)^' (w) are calculated using Eq. 
with Eij>o -» Eij>o and x v( u ) -» X\f{jp,(j) or 
X\p(w) from Eqs. (|A"7)) - (fA"8|) respectively, while A (p) 
is the zero-mode contribution defined by 



A (p) = - [C^f*^^ + cffi*W^(iV< + 1)' 



f 



[qg a) *uipJVi + qg 1} %(JVi + i) 



(A12) 



A R (p) is the UV-renormalization (if any), obtained from 
the last line of Eq. {KQ via A p H = aR (p)M w )- 



1. Special cases 



In this paper, we are interested in the quantities 
AEpp(u;) and Ap'q (w), defined respectively as the 
condensate and perturbation parts of the particular pro- 
jection of Sn and 5m R given in Eq. (A28) of Writing 
AEpp\uj) in the form of Eq. (|A11() we have 

AE^(uj) = A4+>H + AE (p) + AE R (p). (A13) 



Comparing with Eq. (|A4I) . we find that for AE PP ~\u)), 
AE (p) and AP^q (u>), the C P y coefficients are given at 
finite temperature (f(Ni) = 2Ni) by 



C 



(A) 



Y 



(A) 



r (Bl) = y(Bl) r (B2) _ y(B2) 
pij pij ' pij pij 



where the Y„ 



(A14) 

Ypij coefficients are the same as those in 
Eq. ljA7j) and are defined in Eq. (106) of °3]. At zero 
temperature (f(Ni) = 1) these coefficients are given in- 
stead by 



C (A) = Y (A) 



Pi] 



r (Bi) = Y 

pij pi'J 



pi] 
(SI) 



dA 



pi] , 



C 



(B2) 
pij 



Y, 



(B2) 



dBl 

dB2 r 



pij, 



(A15a) 
(A15b) 
(A15c) 



where the quantities dA p ij etc are the changes to the 
integrals A pi j etc defined in Eqs. (107)-(109) of aris- 
ing from the asymmetry in the zero-temperature part of 
Sh(r,uj). They are given by 



dApij = 2V AWo Jdr (u p $q + v p ^ )(viUj - UiVj), 

(Al6a) 

dBlpij = 2^/NqIIo Jdr (u p $q + v p $o)(viUj - UiVj)* , 

(A16b) 

dB2 pij = 2^/N^Uq Jdr (u p % + v p <P )(v*v J - U*Uj). 

(A16c) 



The UV-renormalization AE R (p) is given by Eq. (105) in 
3] and involves the quantity AUq of Eq. I|6|l. For numeri- 
cal consistency, we calculate AUq using an integration up 
to E cut and only take the integral to infinity if we also 
include a semi-classical approximation for higher energy 
states. 

To calculate the condensate number fluctuations 
5Nq(u>) = — Jdr 5h(r,uj) (which are only non-zero for 
modes with m p = p p = 0), we have at finite temperature 



2Ni 



C 



(A) 



C 



(Bl)* 



'N Q Ji 



i], 



C jHj 2 - -V^ohj, 



(A17) 



where Zy and Jy are defined in Eqs. (HO)-(lll) of 
At zero temperature (f(Ni) = 1), the coefficients C, 



and C„; A are unchanged, while C pi ^ becomes 



(A) 

pij 



pij 



c. 



(B2) 



= -2jN n V,, 



o v ij , 



(A18) 



with Vij as in Eq. (100) of |3j. In this case there is no 
zero-mode contribution because the condensate is orthog- 
onal to states of finite energy, and there is also no UV- 
renormalization. 



2. Condensate density fluctuations 

Experiments generally measure moments of the con- 
densate density fluctuations given in Eq. (|16fl . If only 
the single mode 'p' is excited these fluctuations have 
the spatial form u p (r)<I>g(r) + u p (r)<I>o(r). Since the den- 
sity fluctuations are real in the time-domain, we calcu- 
late their projection onto the real part of this function, 
f p {r) = Re[M*(r)$ (r) + v*(t)$>$(t)]. Defining the inte- 
gral 



dr/ p (r)[ Up (r)$*(r) + v p (t)9 {t)], (A19) 



we therefore calculate the quantities 



0**'M or^M 



Nod p P p o(uj) 



dr f p (r)Sn c (r, uj) 



Re(c p ) 
N U Q dp 



SN (uj) , 2b No (w) 



N P pQ (uj) P p o(w) 



(A20) 
Ppo (w) ' 



where c p is defined by the equilibrium limit of Eq. JSJl, 
6at (ijj) is the coefficient defined in Eq. (81) of Ref. |3| and 
we obtain either Q p nc (uj) or lZ s p nc (uj) on the left-hand-side 
depending on whether or not we exclude direct driving 
of the non-condensate in evaluating the expressions on 
the right-hand-side. If either m p ^ or p p ^ 0, then 
SN (uj) = b No {uj) = and G s p n °{uj) and n 5 p n -(uj) reduce 
to the resolvents Q p (uj) and 1Z p (uj) defined in the main 
text in Eqs. (|23|l and (|21|) respectively. We see therefore 
that these quantities are directly proportional to the con- 
densate density fluctuations measured in experiments. 
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3. Non-condensate density fluctuations 

To calculate the relative phase of condensate-non- 
condensate oscillations, or simply to compare the relative 
sizes of condensate and non-condensate fluctuations, we 
also need to calculate the projection of Sh(r, u) onto the 
function f p (r) defined above. For direct comparison with 
the resolvents Q p n "(uu) or 7Z p nc (uj), we actually calculate 



or^fH 



dr f p (r)5h(r,uj) 



N d p P p o(uj) 

where we obtain either G^(ui) or nf {uj) if we exclude 
or include direct driving of the non-condensate respec- 
tively. Thus the total density fluctuation projections are 
proportional to G 5 p n "(uj) + Qf-iuj) or K p n -{uj) + K s p "(uj). 
The phase of the non-condensate oscillations relative to 
those of the condensate can be found as a function of 
frequency from the argument of the complex quantities 
[S£ a (w)]* x e**«(w) or [TLf(uj)]* x TZ s p n -(uj). 

The projection Jdr f p (r)Sh(r, ui) can be found from 
Eq. (|A2|) and has the generic form given in Eq. (|A4|) . For 
numerical convenience we actually calculate the quantity 
2U a Jdr [w*$o + v*$>o}Sh(r, oj), for which the various C pi j 
coefficients are given at finite temperature (/(iVj) = 2Ni) 
by 



C 



(SI) _ 



C 



Cptj = 2\/%Uo J dr {u p $* + Vp^oXmvj + v iUj ), 

(A22a) 

2-y/Wo Jdr 

(A22b) 



(S2 

pij 



(A22c) 

For the zero temperature contribution (f(Ni) = 1) these 
coefficients are modified according to 



a 



(A) 



c, 



pi j pi>3 
(SI) ^ 



r (j4) A- d A ■ ■ 



Cpij +dBl pi j, 



pij Pij PV i 



(A23a) 
(A23b) 
(A23c) 



with dApij etc as in Eq. IJA16I) . No UV renormalization 
is required in this case so A^ n (p) =0. 



APPENDIX B: SEMI-CLASSICAL 
APPROXIMATION 



very accurate even for quite low values of the cutoff. Un- 
fortunately, the situation is much more complicated for 
the dynamic terms, however, and additional approxima- 
tions are required to reduce the expressions to a manage- 
able form • This limits the practical accuracy of the 
method in this case and we have found numerically that 
the approximation we implement for the dynamic terms 
only becomes adequate at quite high energies and must 
be used in conjunction with a large numerical basis (see 
Sec. HHH and Fig. |TJ. 

The semi-classical approximation is described in detail 
m Refs. [UliHi^. The essence of the method is that the 
quantum numbers 'p' labelling a particular quasiparticlc 
wavefunction u p (r) become a momentum label p so that 
the new wavefunction is defined in a single-particle phase 
space (p, r) by 



u p (r) -> u(p,r) = u(p,r)- 



Jp.r/h 



v p {r) -> w(p,r) = u(p,r)- 



V 

Jp.r/h 



V 



(Bl) 
(B2) 



where V is some suitable volume. The exponential factor 
contains the 'fast' dependence of the functions on the co- 
ordinates, while u(p,r) and w(p, r) only vary on a scale 
set by the size of the condensate and the chemical poten- 
tial. The semi-classical solution to the BdG equations is 
then given by 

-2/ \ i i -2/ s <^(P,r) + e(p,r) 
u (p,r) = 1 + v (p,r) = — , (B3) 



u(p,r)w(p,r) 



2e(p,r) 
n {r)U 
'2e(p,r)' 



(B4) 



where no(r) = No\&o(r)\ 2 is the local condensate density 
and e(p, r) and e HF (p, r) are the local quasiparticle and 
Hartree-Fock energies respectively, defined by 

e 2 (p,r) = 4,(p,r) - (n (r)C/ ) 2 , (B5) 
P 2 

e HF (p, r = + y trap (r) - A + 2n (r )U . (B6) 
2m 

Summations over quasiparticle states are replaced by in- 
tegrations over momenta via 



V ^ 



dh 



{2irhf 



(B7) 



We therefore obtain the following semi-classical ap- 
proximations to the equilibrium non-condensate mean- 
fields n(r) and ?ri R (r) needed in the static terms 



n sc {r) 



d 3 p 

(2irh) : 



[(u 2 +v 2 )N{e) + v 2 ]<d(e- E cut ), (B8) 



The calculation of both the static and dynamic terms 
may be supplemented with a semi-classical approxima- 
tion for the quasiparticle modes above the numerical cut- 
off energy to improve convergence. This procedure is 
straightforward to implement for the static terms, and is 



1 



d 3 p n U [2N(e) + l] 

{2-Khf 

d 3 p 



2c 
n U Q 



(2irh) 3 2(p 2 /2m)' 



e(e-E cut ) 



(B9) 
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The last line of the expression for m R (r) contains the full 
UV renormalization of the anomalous average (necessary 
for convergence at zero temperature) and is integrated 
over all energies, not just those above E cut . In practice, 
we divide this integral into two regions, corresponding to 
energies above and below the cutoff. The contribution 
from energies below the cutoff is included in the exact 
basis calculation of m R (r) so that this is well-behaved 
as the cutoff increases, while the part above the cutoff 
is calculated along with the remaining semi-classical cor- 
rection in Eq. (|B9|> . 

The integrals over the direction of the momentum in 
Eqs. (|B8|I and l|B9|) arc trivial so the semi-classical ap- 
proximation for the static terms reduces to the calcu- 
lation of onc-dimcnsional integrals over energy for each 
spatial grid point. In this case, the approximation rapidly 
becomes highly accurate, so if only static terms are re- 
quired in the theory a small numerical basis is sufficient. 

The problem of how to apply the semi-classical ap- 
proximation to dynamic terms of the form of Eq. I|27|) 
has been studied in detail by Giorgini |17| . In a tour-de- 
force calculation, he combined the semi-classical treat- 
ment of the excitations with a Thomas- Fermi approxima- 
tion for the condensate (necessary for consistency at low 
energy 0|) to obtain analytical predictions for energy 
shifts and decay rates for trapped gases in the thermo- 
dynamic limit. Unfortunately the method is limited to 
modes with constant Laplacian because of the nature of 
additional approximations which have to be made to deal 
with the Landau terms. Although these modes include 
those studied to date, we require a method which can be 
used for finite systems outside the Thomas- Fermi regime. 
This can be obtained by adapting the method described 
by Giorgini to the current problem. 

We start by writing the semi-classical approximation 
to the dynamic terms in the form (we focus on AE p D ^ in 
the following) 

&EW(p) = — JJJ ^Wsd 3 Pi d 3 qe- iq s 

9p( r ) K i 3 ( r > r + s, u)g P (r + s) 

(BIO) 

where q = p, — pj , s = r' — r and g p (r) is a smooth func- 
tion of r dependent only on $o(r), u p (r) and v p (r) (these 
are the full numerical solutions and do not involve any 
semi-classical approximation). The kernel Kij(r, r + s, oj) 
contains all the slow dependence of the intermediate 
states, i.e. the energy denominators, population factors 
and the slowly-varying part of the quasiparticle wave 
functions. The exponential factor e~ iqs contains the 
rapidly varying part of the quasiparticle wave functions 
and acts like a delta function in both position and mo- 



mentum, so that an expansion in powers of s is appro- 
priate. The leading order corresponds to setting s = in 
the kernel and integration over s in the exponential then 
leaves a delta function of q. The double integrals over 
position and momentum therefore collapse into single in- 
tegrals as for the static shifts, and we obtain 

AE s?Hp) = T^U / f<?v<P Vl |. 9p (r)| 2 ^(r,r,c). 

(Bll) 

The main difficulty with applying this approach is the 
pole in the kernel for the Landau terms where ~ 
(Nj — Ni)/(uj — €j, + €j). This invalidates the assump- 
tion that the kernel is slowly varying and restricts the 
analytical calculation to modes with constant Laplacian 
|17|. However, we are only interested in the semi-classical 
approximation above a large cutoff energy, and we there- 
fore expect that the result will be approximately inde- 
pendent of frequency in the range of interest. We there- 
fore set u> — in the kernels, which allows us to replace 
the badly-behaved factor (Nj — Ni)/(uj — ti + Cj) with 
the derivative —(dN/dE)\ Ci . The final form of the semi- 
classical approximation that we actually use for AE p D ^ 
is therefore 

(l + 2iV)([y^ ) (r)] 2 + [Y p ( f ) (r)] 2 ) 

4e(p,r) *> 
(B12) 

where the coefficients Ypn(r) are the integrands of the 
coefficients Y™,- appearing in Eq. I|A7|I and defined in 
Eq. (106) of y| but with the exact quasiparticle wave- 
functions replaced by their slowly varying counterparts, 
i.e. Ui(r) —> u(pi,r) etc. We have found numerically 
that the approximation of Eq. (|B12(I is much less accu- 
rate than the semi-classical results for the static terms 
but becomes adequate when a high cutoff energy is used. 
The effect of the semi-classical approximation on our nu- 
merical convergence can be seen in Fig. ^ 

We have not used a semi-classical approximation for 
the dynamic term AP^'(w) because the long length 
scale of the external potential invalidates the approach. 
This is not particularly serious, however, as there is no 
delicate cancellation of terms in the numerator of the 
response function as there is in the denominator so the 
accuracy requirement is greatly reduced. In addition, the 
use of the asymmetric summation technique described in 
Sec. IIII El removes any need for a semi-classical approxi- 
mation for this term. 
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